function f = z_high1(zh,beta,k,R,q_bar)
% Upper bound of branches leading to bunching: fixed
f = exp(zh + beta*R - 1)/beta - k -(R - (q_bar - zh)/beta)*exp(q_bar); 
end
